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Abstract. The stochastic motion of a two-dimensional vesicle in linear shear flow is studied at finite 
temperature. In the limit of small deformations from a circle, Langevin-type equations of motion are 
derived, which are highly nonlinear due to the constraint of constant perimeter length. These equations 
are solved in the low temperature limit and using a mean field approach, in which the length constraint is 
satisfied only on average. The constraint imposes non-trivial correlations between the lowest deformation 
modes at low temperature. We also simulate a vesicle in a hydrodynamic solvent by using the multi-particle 
collision dynamics technique, both in the quasi-circular regime and for larger deformations, and compare 
the stationary deformation correlation functions and the time autocorrelation functions with theoretical 
predictions. Good agreement between theory and simulations is obtained. 

PACS. 87.16.Dg Membranes, bilayers, and vesicles - 87.15.Ya Fluctuations - 67.40. Hf Hydrodynamics in 
specific geometries, flow in narrow channels 



1 Introduction 

The dynamics of soft objects such as drops, capsules and 
cells in flow represents a long-standing problem in sci- 
ence and engineering, but has received increasing interest 
recently, in particular due to its relevance to biological, 
medicinal and microfluidic applications. This problem is 
challenging from a theoretical point of view, because the 
shape of these objects is not given a priori, but determined 
dynamically from a balance of intcrfacial forces with fluid 
stresses. Improved experimental methods have revealed in- 
triguing new dynamical shape transitions due to the pres- 
ence of shear flow. The phenomenology of the dynamical 
behavior depends distinctively on the specific soft object 
immersed in the flow with fluid bilayer vesicles and elastic 
microcapsules as the most prominent classes. 

Fluid bilayer vesicles assume a stationary tank-treading 
shape in linear shear flow, if there is no viscosity contrast 
between interior and exterior fluid 1]. If the interior fluid 
or the membrane becomes more viscous, a transition to a 
tumbling state can occur @, 1, i, i, i, 0]. Tank-treading 
was observed experimentally in infinite shear flow 0, [j| 
and for vesicles interac ting with a rigid wall [13, EH, where 
a dynamic lift occurs 0, O E3, O • The tank-treading 
to tumbling transition was observed for the first time con- 
vincingly in an experiment only very recently [T^ |. In ad- 
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dition to the tank-treading to tumbling transition, an os- 
cillating motion was predicted theoretically [17 1 a nd ob- 
served experimentally [l6j and in simulations 18|. This 
type of motion has alternatively been called vacillating- 
breathing [ijj , swinging [l8| , or trembling [lj| . The 
theoretical description has been extended recently beyond 
first order in the shear rate [l8|, QjJ, H(| • 

At finite temperature, stochastic fluctuations of the 
membrane due to thermal motion affect the motion of the 
object. Due to the dissipative nature of the hydrodynamic 
interactions, vesicles in shear flow form a non-trivial model 
system for studying non-equilibrium stochastic dynamics. 
Since the effect of thermal noise on the transitions between 
the different modes of motion in general is a challenging 
task, in this paper we concentrate on the stochastic mo- 
tion in the stationary tank-treading state. Our theoretical 
approach is similar to that of Ref. (UJ, where stochastic 
equations of motion were derived for quasi-spherical vesi- 
cles. 

Most numerical methods solving the equation of mo- 
tion of vesicles or capsules [H, [|| operate in the absence 
of thermal forces. An exception, which naturally includes 
thermal noise, is multi-particle collision dynamics (MPC), 
also known as stochastic rotation dynamics (SRD) [22|, 
HI, H3, Hl|- In this method, the fluid part is modeled 
on a particle rather than a continuum level. The micro- 
scopic equations of motion for the effective fluid are cho- 
sen to be evaluated efficiently on the one hand, and on 
the other hand to lead to the correct macroscopic hydro- 
dynamics. This method has successfully bee n applie d to 
flow around rigid objects [2f| [27j , polymers [28l . l2fjj and 
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viscous vesicles @, @, G]| . We employ the MPC simula- 
tion method to compare our theoretical predictions of cor- 
relation functions, inclination angles, and tank-treading 
frequencies with simulation data of vesicles. In order to 
obtain good statistics, we focus here on two-dimensional 
(2d) vesicles with a linear boundary. 

The paper is organized as follows: After formulating 
the problem in section [21 we develop nonlinear stochastic 
equations of motion for quasi-circular vesicles in section 
[3l These are solved approximately using a mean field ap- 
proach and a low temperature expansion in section 0J We 
also present the 2d version of the deterministic Keller- 
Skalak theory [3l| in section O which takes into account 
the influence of the vesicle shape on the flow. The simu- 
lation method used is discussed in section [6] Finally we 
compare the calculations with simulation data in section 
[7] and discuss our results. 



2 Problem formulation 

We consider a model 2d vesicle immersed in a fluid of 
viscosity ?7 ou t with a Id membrane boundary surrounding 
a fluid of viscosity 7/; n and at finite temperature T. Due to 
the incompressibility of the membrane and of the enclosed 
fluid, the area Aq and the length Lq of the membrane 
are constants. The membrane resists deformation with a 
bending rigidity n, which is defined rigorously below in 
section 12.11 The fixed area defines a length scale 



(1) 



which can be used to define a number of dimcnsionlcss 
quantities. In the following, we use the excess length 



J»0 



2tt 



and the dimensionless viscosity contrast 

'/in 



A = 



?7out 



(2) 



(3) 



Alternatively one can derive a length R* = Lq/(2tt) from 
the length constraint, and use it to define a reduced area 
A* ee A /(nR* 2 ). The reduced area is connected to the 
excess length by 



.4* 



1 



A 
2^ 



(4) 



In a quiescent fluid, thermal stochastic forces acting on 
the membrane lead to a fluctuating shape, where the prob- 
ability of any specific deformation can be calculated using 
the Boltzmann weight corresponding to the deformation 
energy H[r]. If an external flow field v°° is switched on, 
the system ceases to be in equilibrium, and the statistical 
weight of a deformation cannot be calculated a priory us- 
ing Boltzmann weights. We first derive the force balance 
governing the motion of a vesicle in stochastic Stokes flow, 
before we simplify the equations of motion in the limit of 
small deformations from a circular shape. 



2.1 Constitutive equation of the membrane 

We employ conventions of differential geometry following 
Ref. [32|. The shape of the vesicle is given by the shape 
function r(s), where < s < L denotes the arc length. 
The tangent vector t(s) ee dr(s)/ds is of unit length. The 
unit normal vector n(s) is defined to point to the outside 
of the vesicle, and the orientation is chosen such that the 
pair (n, t) forms a right handed system. The curvature 
k(s) is defined via the relation 



dt(s)/ds = -fc(s)n(s). 



(5) 



The 2d analog of the bending energy of a certain mem- 
brane deformation is given by the Helfrich term [33j 



n K [r] - - 2 



&sk(sf 



(6) 



which corresponds also to the bending energy of a semi- 
flexible polymer [34]]. Note that for 3d vesicles, a sponta- 
neous curvature Co can appear in the bending energy for 
intrinsically asymmetric monolayers or asymmetric liquid 
environments. In 2d vesicles, we can ignore the sponta- 
neous curvature, since it shifts the bending energy only 
by a topological constant, much like the Gaussian curva- 
ture contribution to the curvature energy can be ignored 
in 3d. 

All deformations of the vesicle must preserve the length 
L. In addition, the fluid membrane is locally incompress- 
ible. This is ensured by introducing the tension a{s) as a 
Lagrange parameter. The total energy thus reads 



H[r\=H K [v]+ [ dsa{s). 
Jo 



(7) 



From the Euler-Lagrange equations we can deduce the 
force acting on the membrane 



t (a' + 2nkk') + n Qfc 3 - nk" - ka 



! = SH[r] 
Sr 

Here the prime denotes a derivative with respect to the 
arc length s. 



2.2 Stochastic Stokes flow 

The elastic forces given by Eq. ([8]) are balanced by hydro- 
dynamic forces mediated by the surrounding fluid. The 
motion of the fluid and the vesicle is only considered on 
time scales on which the fluid is incompressible, i.e. 



0. 



(9) 



The length and time scales in typical experiments and sim- 
ulations is such that the Reynolds number is very small. 
We only consider fluctuations on time scales on which the 
incrtial term in the Navier-Stokes equation can be ne- 
glected. The velocity field v of the fluid is then subject 
to the steady stochastic Stokes equation [35| 



Vp + r] a Av + V -s = 0, 



(10) 
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where the thermal stress tensor s(x) is assumed to be a 
Gaussian random variable with zero mean and correlations 

(s ik ) = 

(s. lfc (xi, ti)sz m (x 2 , t 2 )) = 2k B Ti] a S(-K 1 - x 2 )<5(£i - t 2 ) 

Ohm 

(11) 

Here a € {in, out} indicates the inner or outer fluid. In- 
stead of calculating the stochastic velocity field of the flow, 
we only calculate the deterministic part of v: 



Vp + r) a Av = 0. 



(12) 



At the vesicle membrane we must have force balance be- 
tween the deterministic and stochastic part of the hydro- 
dynamic force and the elastic forces 



[T-n]l 



s • n = 0. 



(13) 



Here, T denotes the deterministic hydrodynamic stress 
tensor with Cartesian components 



T ik = ~pS ik + Va[diV k + d k Vi 



(14) 



Far away from the vesicle the velocity field assumes the 
externally given values 



v(x)^v°°(x), |x|-»oo, 



(15) 



which is ensured by separating an induced part from the 
velocity field 



v = v°°+v ind , (16) 



and requiring that the induced part drops to zero far away 
from the vesicle. Assuming no-slip boundary conditions, 
the vesicle is advected by the flow, which implies 



d t r(s,t) =v(r(i,i),t). 



(17) 



Here the dynamics still depends implicitly on <r(s), which 
has to be chosen such that s remains the arc length, ensur- 
ing incompressibility. Eqs. (|8ll2ffl7|) determine the stochas- 
tic motion of the vesicle. 



3 Quasi-circular approximation 

These equations can be simplified considerably if we re- 
strict ourselves to vesicle shapes close to the circle. We 
parameterize the shapes as a function of the polar angle 



r(^) = iJ o e r (^)(l+«(0)), 



(18) 



and consider small distortions u. The deformation ampli- 
tude u{4>) is a real periodic function of (f> and can therefore 
be expanded into complex Fourier modes 



u{4>) 



m— — oo 



exp(zm^)) 



(19) 



For comparison with simulation data described below, an 
expansion into a real Fourier series is advantageous. We 
therefore also employ the expansion 

oo oo 

u{4>) = ciq + a m cos(m(f)) + b m sin(m</>). (20) 



m— 1 



m— 1 



The real Fourier coefficients a m , b m are connected with 
the complex Fourier coefficients u m via (to =/= 0) 



2 \P"m ^m) - 



(21) 



Area conservation fixes uq in terms of the other u m 

tin = ) >> m \ 2 . (22) 



?b E K 



This relation will be used throughout the paper, and from 
now on sums over to exclude the to = term. The contour 
length L of the membrane is calculated to second order in 
u to be 



L = 2nR + ^ ^(TO+l)(TO-l)| Um | 2 



Hence the excess length A reads 



TO 2 - l) \u m \ 2 . 



Finally, the local curvature k evaluates to 

R Q k(<f>) = l-u"{<t>) = l+Y J m 2 u m - Mlm ^ 

m/0 



2tt 



(23) 



(24) 



(25) 



This leads to the bending energy (ignoring constant terms) 



= iii E K - 1) K - 3 / 2 ) 



(26) 



We now add the global length constraint (|24|) with a La- 
grangian multiplier 



no/R\ 



(27) 



to the quadratic part of the bending energy. This leads to 
a quadratic expression for the total energy 



(28) 



with 



E m {a) = (to + l)(m - 1)[to 2 - 3/2 + a}. (29) 



The bending forces © are determined by the deformation 
amplitudes u m and by the instantaneous tension 



= £ ( ° + E ct ' 



exp(«TO0) 



(30) 



The homogeneous tension na/Ro has already been in- 
cluded into the energy (|28|) . 
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3.1 Velocity field 

In polar coordinates, the general solution of Stokes' equa- 
tion can be expanded into the fundamental modes [36| 



m 



V( r ±M exp(zm<£))/-\ 
0, 



p -± = 1 _ (M±H r 2 V(r±l" l l exp(zm0)) 



Prrl = V<*r =t| 



2(l±|m|)v / 27r V 2 l m l 

rr ±M exp(im0)) 



(31) 



2tt. 



In this representation the cases \m\ — 1 and \m\ = 
are special and have to be treated separately. They corre- 
spond to constant flow and rotational flow, respectively. 
The deeper reason why these are special cases is the Stokes 
paradox [36[ . It follows from the boundary conditions that 
the induced velocity field on the inside must be composed 
of "+" modes, and of "— " modes on the outside. 

The corresponding hydrodynamic stress tensor reads 
in (r, (f>) components 



p<2> ± 



2r\r 



±|m|-2 ex PO^ 



2tt 

|m|(|m| =p 1) im(±|m| — 1) 
im(±|m| — 1) — |rn|(|m| =p 1) 



(32) 



and 



-,p.± 



rjr 



±|m| exp(im<? 



2V2tt 



±|to| - 
im 



2 im 

tM - 



(33) 



We can now express the 2d Oseen tensor in spectral com- 
ponents. The radial and polar components of the fluid 
velocity and hydrodynamic force at the reference circle 
are expanded into Fourier modes analogous to the expan- 
sion (fl9|) . The velocity field at the reference sphere to- 
gether with the boundary conditions uniquely determines 
the expansion (f3Tj) . From the spatial velocity field the hy- 
drodynamic force f = T ± ■ n can be calculated, leading 
to 



jr,ind 
J m 

j*0,ind 
J m 



Vin + Vont 

Ro 



2\m\ 2i sign(m) 
-2% sign(m) 2\m\ 



,r,ind 

m 

m 

(34) 



there is a jump in the traction 
exp(zm 



r,oo 
(b.oc 



E 



2tt 

2|to|(|to| - 1 



(r)in - Vout) 



2im(\m\ - I) I m 



m\/2 - 1 
im/2 



poo 
m 

(36) 



For the specific case of external linear shear flow 



iye x = { A f/2){ye x + xe y ) - (j/2)(xe y - ye x 



-(ijV2n/8) 



v, — v 



- (7/2)re , 



we can read off the only non- vanishing components 



iV27T7/8, 



Q = -7/2. 

We will also use the dimensionless shear rate 



and vorticity 



X = 7 



n = n VoutR3 ° 



x 
2' 



(37) 



(38) 



(39) 



(40) 



3.3 Incompressibility condition 

The flow at the vesicle membrane is subject to the incom- 
pressibility condition D t ^/g — 0, which can be cast in the 
equivalent form t • 9^v(r(0)) = 0. To leading order in the 
deformation, this condition reads 



v r (R ) + d<j,v<j,(Ro) = 0. 



(41) 



Separating the induced flow from the external flow, we 
have in Fourier components 



3.2 External flow 

In the absence of the vesicle the applied external flow must 
be regular everywhere. Therefore apart from constant flow 
and constant rotation only the "+" modes contribute in 
the expansion (|3ip . To avoid the intricacies of the Stokes 
paradox, we neglect the possibility of constant flow. A 
general expansion of the external flow therefore reads 



v *.- 

m m 



(35) 



The last term in this expansion corresponds to rotational 
flow with the vorticity ft. For a finite viscosity contrast 



0,ind 



= M(H-i)4 ro| - 1 as 



4(|m| + 1) 



(42) 



Using this relation, we can eliminate and obtain 
/ m = 2— (A + l)^— , 



r.ind 



+ #~2(|m| - 1)^ [|m|(A - 1) + (A + 1)] + P= 



(43) 
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3.4 Equation of motion 

Neglecting the thermal fluctuating forces for the moment, 
we can derive a deterministic equation of motion. The 
force balance leads to 



u 

with 
and 



ind,r 



-(K/r)outRl)r m E m (a)u m + B m $^R Q 



771 



2(A+ l)(m 2 



B, 



777 2 (A-1) + H(A+1) 

(A + 1)(H + 1) ■ 



(44) 
(45) 

(46) 



From the induced velocity, we obtain the radial component 
of the full velocity field 



„r,ind 



r.ind 



(47) 



, , . «m *™ + C — v m"~ + \ m Wm- 

The advection equation then reads (cf. Ref. [2l|) 

d t u m = iflmu m - {K./r/ out Rl)r m E m (a)u m + D m <P% (48) 
with 



D„ 



2m 
A+l 



(49) 



At non-zero temperature, thermal forces must be taken 
into account in the force balance. The deterministic equa- 
tion of motion (|48[) then becomes a Langevin equation 

d t u m = if2mu m - (K/7] out Rl)r m E m (a)u m + D m <P^ + ( m . 

(50) 

The form of the thermal noise £ m can be obtained directly 
from the noise term in Eq. (|13p . It is much easier, however, 
to determine £ m from the Einstein relation, which must 
be valid in equilibrium. We assume that the equilibrium 
noise is valid also for non- vanishing shear flow and choose 

(Cm(*)Cm'(*')) = 2(k B Tr m / Vout R 3 )S rn, — m' 

5{t-t'). (51) 

Eq. ([50)) is the correct stochastic equation of motion for 
the vesicle deformation modes u m . The tension a is at 
each instance determined such that the length constraint 
(HU) is fulfilled. Taking the time derivative of Eq. (gj) and 
using Eq. (|50p. we can solve for the tension 



m^O 
K 

?7outi?o 



(ttt 2 



iflm\u n 



r m ^; m (o)|u m | 2 + D m u* m <p^ + u* m c, 



(52) 



When this expression is inserted back into Eq. (f50|) . the 
resulting noise term becomes dependent on the instanta- 
neous values of the u m . While such non-linear noise terms 
hold interesting physics, we first concentrate on tractable 
approximate solutions to the stochastic equation of mo- 
tion. 



4 Approximate solutions 

4.1 Mean-field treatment 

At finite temperature, higher-order modes are excited by 
stochastic thermal forces and therefore cannot be neglected. 
The full non-linear set of Langevin equations (|50|) in com- 
bination with the expression (|52[) for a is too complex to 
admit a general solution. We can, however, gain further 
insight in the tank-treading regime using a mean-field de- 
scription. We replace the fluctuating tension a in Eq. (|50l) 
by a constant, which has to be determined self-consistently 
from the length constraint. The Langevin equations (|5"0")) 
then become linear and decouple. In the stationary state, 
only the m — 2 deformations have a finite mean, 



(u 2 = 



D 2 



r 2 E 2 {(j) + ix 



(53) 



On average, the vesicle is elliptical. As a measure of the 
deformation from the circle we define the Taylor deforma- 
tion parameter 

^Irl- < 54 > 

where L and S denote the long and short axis of the ellipse. 
In the mean-field treatment we have 



D 



[(5/2 + (t) 2 + 9x 2 (1 + A): 



,i/2- 



The inclination angle is obtained from Eq. ([53 


The deviations from the mean 



1 5/2 + 0- 
- arctan — ; — . 

2 X(l + A) 



(55) 



(56) 



5u m = u m - (u m ) (57) 

obey the homogeneous Langevin equation 

d t Su m = imf25u m - (n/r] out Rl)r m E m (cr)6u m + ( m . (58) 

The stationary noise correlations are best evaluated using 
a time Fourier transform 



5u{uj) 



dt e~xp{—ibjt)8u{t), 



(59) 



leading to 

iu)5u m = -imf26u m - (n/r] out Rf,)r m E m (a)6u m + (, 
We can solve for 5u m (uj) and obtain the correlations 

2k B Tr m /r/ out R% 



(60) 



(5u m (u)5u- m (-u))) 



(lj + mSI) 2 + 



Kr m E,„ 



(61) 
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We have left the cr-dependence of E m implicit for clarity. 
The time correlation function becomes {At > 0) 

/d 
— exp(iivAt)(8u m (<jo)5u- m (-ui)) 



fluctuating parts to the excess length can be determined 
analytically to be 



k B TR Q 

K,E m 

x exp 



( Kr m Ern 

V W«0 



Q] At 



with the stationary equal-time correlations 

k B TR a 



(5u m (t)5u- m (t)) 



nE m {a) ' 



(62) 



(63) 



The amplitudes u m with different m are uncorrelated at 
all times. Comparison with simulation data is easier us- 
ing the real Fourier coefficients (|2"0|) . The corresponding 
correlation functions read 

(Sa m (0)Sa m {t)) = (b m {0)b m (t)) 
k B TR 



■ exp 



KT m E n 



lTKE m ""^ V VautR'o 

(5a m (0)Sb m (t)) =-(6 TO (0)o m (t)) 
k B TR Q 



f-f ) cos(?7J7i/2), 



irnE n 



(-SIT*) Mmjt/2), 
(64) 



and 



(Sa m {0)5a m (0)) = (6b m (0)6b m (Q)) = 
(5a m (0)5b m (0)) = 



k B TR 
■KKE m [a) (65) 



The fluctuating u m contribute to the excess length accord- 
ing to Eq. ([24]) . Although the length constraint cannot be 
obeyed exactly with a constant tension, we determine a 
such that the constraint is fulfilled on average. The 
total excess length has a systematic and a fluctuating part 



A = A{a) + A m(cr), 



(66) 



with 



A(a) 



E(- 2 

m>0 
37T 



n 2 \<p°°\ 2 

ml m I 



x 2 



i 2 f2 2 



(67) 



2 (5/2 + a) 2 + x 2 (l + A) 



2 ' 



and 



A m {a) = (m 2 - l)(\6u m \ 2 ) 



k B TR 



n{m 2 -3/2 + a) 



(68) 



Thus cr is determined implicitly by the solution of Eq. (|66[) . 
For future reference, we note that the contribution of the 



E, i \ k B TR 
../ m(CT)= Wf(^ 



v^(7 - 6(j) 



(4cr 2 -8a + 3) 

+7rv / 3 — 2^(2cr - f ) cot - cr) 

(69) 

While this expression is exact, its behavior as a function 
of a is not obvious (for example, the "singularities" at 
(7 = 1/2 and cr = 3/2 are only apparent). We therefore 
give the leading asymptotic behavior 



^ , fc s Ti? fl/((7 + 5/2) + 25/48 ct 
> An(cr) « < 1/2 

m>2 K v ; 



-5/2 

00. 

(70) 



4.2 Zero temperature 

At large shear rates, nearly the entire excess length is 
stored in the systematic part A. As a crossover shear rate 
Xcj w e can define the shear rate at which the two contri- 
butions in condition ([66)) become equal 



Afa Xc) 



m>2 



(71) 



This set of equations must be solved numerically for each 
A. In the limit x ^ Xc we can ignore the thermal forces. 
In this case, the equation of motion (f50|) becomes the 
deterministic Eq. (|48|) . and the tension is determined by 
Eq. (53H with Cm = 0. 

We can easily obtain the stationary state from c>tM m = 
0, i.e. 

= Tfout-Rp A" ^oo / 72 x 

l m E m (<T ) + im$2 

The homogeneous tension cto is determined from the length 
constraint (f2~4"|) . In the case of constant linear shear flow, 
only the m = ±2 components are non-zero and are equal 
in magnitude. The length constraint thus reads |u±a| = 
{A/3) 1 ' 2 , or 



A = 3 



Do 



2w X 2 



r?E 2 

3tt 



4f2 2 64 



(73) 



2(1 + A) 2 (5/2 + cr ) 2 /(9(l + A)) 2 + X 2 



The homogeneous tension in the stationary state is thus 
given by 



an = -5/2 + 3 X (l + A) 



3- 



-,1/2 



_2A{1 + A) 2 
Ei vanishes at a critical viscosity ratio 



- 1 



A i 71 
C ~\2A 



1. 



(74) 



(75) 
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This corresponds to a tank-treading to tumbling transi- 
tion, as can be seen when we allow for time-dependent ao~. 
In linear shear flow, only the m — 2 modes are excited. In 
the long time limit we can therefore assume that all other 
modes have decayed. In analogy with the 3d treatment 
[l7| . we can write u 2 in polar form 



where the noise term 



u 2 = (Z\/3) 1/2 exp(-2i6>), 



(76) 



where is the inclination angle of the vesicle with respect 
to the shear direction. Taking the real and imaginary part 
of Eq. (|48|) gives the familiar Jeffery's equation (3?| 



(9 = 7 



/3tt 



2 (A+1)V2Z 



cos(26>) 



(77) 



For A < A c , Eq. (|77|) admits two stationary solutions, of 
which only the positive is linearly stable 



_ 1 ( {\ + l) 2 2A \ 1/2 

°° =2 arCC ° S ( 3n ) ■ 



(78) 



This corresponds to stationary tank-treading motion, where 
the tank-treading frequency at zeroth order is given by the 
external flow 



,QC _ 7 



For A > A c , the right hand side of Eq. ((77)) is always neg- 
ative, and the vesicle starts to tumble. In two dimensions, 
no analogy to a swinging motion (cf. Refs. [13, H, El, 
fl9| ) exists, since the volume and length constraint already 
uniquely determine the shape of an ellipse. 



4.3 First-order correction to the large shear-rate limit 

In the mean-field approach the tension a is assumed con- 
stant and all modes fluctuate independently with ampli- 
tudes given by Eq. ([63]) . In this picture, the length con- 
straint is not fulfilled rigorously but only on average. For 
strictly enforced length constraint the tension must fluc- 
tuate according to Eq. (|52[) . which induces correlations 
between the deformation amplitudes. While this general 
effect is worth considering in its own right, here we con- 
centrate on the much simpler large shear rate (or low tem- 
perature) limit as a perturbation of the deterministic so- 
lution. 

At T = 0, the whole excess length A is stored in the 
|m| = 2 mode. Perturbing the modulus of the amplitude 
\u 2 \ alters the excess length A to first order and is pro- 
hibited by the constraint j2]) . Perturbing the other modes 
alters A only to second order. At low temperature, we can 
therefore assume the polar decomposition (f76|) . Taking the 
real and imaginary part of the equation of motion (|50"|) . 
we arrive at a Langevin equation for the inclination angle 



= 7 



1 



2 (X + 1)V2A 



cos(26>) 



(80) 



Im[C 2 exp(2i0)] 



is Gaussian and delta-correlated 

3 



mm) 



AA 



r 2 k B TS(t). 



(81) 



(82) 



In the stationary regime fluctuates around the mean 
value 

= O + A0, (83) 



where 0q is given by Eq. (|78|) . For small A0 we can ex- 
pand Eq. ([50)) to obtain 



A0 



3AX+1 

This implies the stationary correlations 



sin(2<9 )^<9 + £. 



(84) 



(A0 2 ) = 



RoksT 1 
nxA 1 / 2 8 



y-(A + l)M 



1-1/2 



3fcsTi?o 



8kAE 2 ((to) 
(85) 

where we have used Eq. ([75)1 . For small A we read off 



(79) (A0 



Q2\l/2 



1 



1297T 



1/4 (k X A^\- 1/2 



\Rok B T 



0.20 



ixA 1 / 2 



Rok B T 
(86) 

Finally, we calculate the fluctuations of the Fourier modes 
a 2 , b 2 • The polar expansion (|76[) implies 



a 2 



62 = 




(87) 



We derive the correlation functions of the m = 2 modes 
from Eq. (|85|) to be 

(Sa 2 6a 2 ) = kB J^° , cos 2 (26> ), 

7TKh 2 (a ) 

(Sa 2 6b 2 ) = kB J R ° , cos(20 o ) sin(20 o ), (88) 



(5b 2 Sb 2 



ttkE 2 ((jo) 
k B TR a 

TTKE 2 (a ) 



sin 2 (26> ). 



5 Keller-Skalak theory 

In the theory of Keller and Skalak [3l| , a three-dimensional 
vesicle is assumed to have a fixed ellipsoidal shape 



(xi/ai) 2 + (x 2 ja 2 ) 2 + {x 3 /a 3 y 



1. 



(89) 



where the ai are the semi-axes of the ellipsoid, and the co- 
ordinate axes Xi point along its principal directions. The 
x\ and x 2 axes, with a\ > a 2 , are chosen to lie in the xy 
plane and are rotated through an angle with respect 



-1/2 
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to the x and y axes. The components of the undisturbed 
shear flow are (jy, 0,0). The velocity field at the mem- 
brane is assumed to be 



v = w tt (-(ai/a 2 )x 2 , (a 2 /ai)xi,0) , 



(90) 



where u>^ is a parameter having the dimensions of a fre- 
quency. The energy supplied by the external flow has to 
be balanced with the energy dissipated inside the vesicle. 
The motion of the vesicle derived from this energy balance 
reads l3fl| 



with 



B = 



f = -j + £cos(2 ), 



7 r(l-r 2 2 ) 2 [z 2 (l-A)-2]-8i 



(91) 



7 f U ~ rj 
1 + r.U 2(1 



(l-r§)Ml-A)-2] 



} (92) 



and 



,KS 



r 2 (l + r|) 



(l-r 2 ) 2 [z 2 (l-A)-2]-8r: 
The factors appearing in Eqs. (|9"Tj) - (l9"3"j) are given by 



(93) 



r 2 = aij a\ 



and 



_ -1/3 -1/3 



r 3 = a 3 /ai, z 2 = g 3 (a 2 + a 2 .), 



Q'2 



_ 2/3 -1/3 



, «3 



_ -1/3 2/3 



(94) 



.93 



(a 2 + s)- 3/2 (a 2 + s)-V 2 (aj + s^ds. (95) 



For B > 7/2, we obtain a steady tank-treading angle 



= — arccos ( — — ) 
2 \2B 



(96) 



in two consecutive steps, streaming and collision. In the 
streaming step, particles move ballistically, 



r i (t + At s ) = r i (t)+v i (t)At s 



(97) 



For the collision step, the system is divided into the cells 
of a regular square lattice of mesh size a. Each of these 
cells is the interaction area where an instantaneous multi- 
particle collision occurs, which changes particles velocities 
as 1 22] 



Vi(t + At s ) = u(t) + JJ[vi(t) - u(i)], 



(98) 



where u is the average velocity of the colliding particles in 
a cell. The velocity field u is considered to be the macro- 
scopic velocity of the fluid and it is assumed to have the 
coordinates of the center of the cell. Q denotes a stochastic 
rotation matrix which rotates, with equal probability, by 
an angle of either +a or —a. The collisions are performed 
simultaneously on all the particles in a cell with the same 
rotation f2, but Q may differ from cell to cell. The local 
momentum and kinetic energy are conserved under this 
dynamics. The kinetic energy of particles fixes the tem- 
perature ksT, where ks is the Boltzmann constant, via 
the equipartition theorem. 

It was shown in Ref. [38| that a proper description of 
hydrodynamics in MPC requires large Schimdt numbers. 
This can be accomplished by choosing a mean-free path 
I = At s y/kBT/m s , which is small compared to the cell 
size a. It is known that a value of I much smaller than a 
breaks the Galilean invariance (39| and that this problem 
can be solved by applying a random sh ift p rocedure (39l |. 
The viscosity of the solvent fluid is [13, Hif 



/ " 
2a . 
a 

'ui 



(n c 



) sin 2 a 



l + e- 



y/m s k B T 



(99) 



We calculate the inclination angle and the tank-trading 
frequency oj§ by adapting the Keller-Skalak theory to 
two dimensions. We numerically solve Eqs. (j9"Tj) - ([9"3")) in 
the limit r% — > +00 keeping r 2 finite, which formally cor- 
responds to an ellipsoid with an infinite semi-axis in the 
z direction. 



6 Simulation method 

A 2d vesicle model system was simulated using the multi- 
particle collision (MPC) dynamics [22j,l5j,l25|. In this method 
the fluid is not treated on a continuum level, but rather 
by a stochastic dynamics of effective fluid particles. 



with particle density p 
tides per cell. 

In order to enforce shear 



lc )(l - cos a) 

J a 

n c m s /a 2 and number n c of par- 



flow, we place our system 



6.1 Solvent dynamics 



of size L x x L y between two horizontal walls. The upper 
and the lower walls slide along the x direction with ve- 
locities v wa u = (v wa u,0) and —v wa u, respectively, with 
Vwaii > 0. Periodic boundary conditions are used along 
the x direction. Along the y direction, we use a mod- 
ified bounce-back boundary condition which consists in 
requiring that particles hitting the walls change their ve- 
locities according to — * 2v wa u — v^. Together with vir- 
tual particles in partly filled cells at walls, this describes 
no-slip boundary conditions very well [2a , |27| . A linear 
flow profile (u Xl u y ) = (jy, 0) is obtained with shear rate 
7 = 2v wa ii/L y , with the walls placed at y = ±L y /2. The 
relative velocities in the collision cells are rescaled after 
each time step At s in order to keep the temperature con- 
stant in the (driven) system. 



We consider a two-dimensional system made of N s identi- 
cal particles of mass m s whose positions r-j(i) and veloc- 
ities Vi(t), i — 1,2, ... , N s , are continuous variables. The 
time is discretized in intervals At s . The evolution occurs 



6.2 Membrane model 

The vesicle membrane is modeled by connecting N p beads 
of mass Trip successively with bonds into a closed ring. 
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Neighboring beads along the closed chain are connected 
to each other with the harmonic potential 



C^bond = "7T 



N p 

£ 

(=1 



(100) 



where fc/j is a spring constant, is the position vector 
of the z-th bead, and tq is the average bond length. The 
bending energy (|26|) is modeled on the discrete level by a 
bending potential 



N p 

Ubend = — VVl ~ COS A), 



(101) 



where Pi is the angle between successive bonds. The fluid 
modeled with the MPC method is compressible. To en- 
force the area constraint in the presence of thermal and 
hydrodynamic forces, we add a constraint potential 



k A {A-A f 



(102) 



6.3 Coupling of membrane and solvent dynamics 

The membrane-solvent interaction must prevent solvent 
particles from crossing the membrane and enforce no-slip 
boundary conditions on the membrane. Therefore we place 
hard disks centered on the membrane beads. The disk ra- 
dius r p is set in order to ensure overlapping of disks and 
a complete coverage of the membrane. The exchange of 
momentum between the solvent particles and the mem- 
brane occurs in the following way. After updating beads 
positions and velocities via molecular dynamics (MD), we 
freely stream all the solvent particles. We then execute 
bounce-back scattering between solvent and membrane 
disks only when a solvent particle j and a disk i satisfy 

the conditions jr* - r,-| < r p and (r; - r^) • (v, - Vj) < 0. 7.1 Stationary deformations 

This means that if the two collision partners i and j over- 
lap and move towards each other, then their velocities are 
updated according to 



the membrane until the inner density increases so that 
an expansion force compensates the compression one. It 
is straightforward to show [6j that the density increase is 
Apj ' p = 2r p /R* where 2r p is the effective membrane thick- 
ness. This requires that Rq is large enough compared to 
the disk radius r p to reduce such compression effects. The 
number of solvent particles placed inside the vesicle fixes 
an average area. However, since the MPC fluid is com- 
pressible, shear and bending rigidity effects may change 
the area A. For this reason the constraint potential (|102ll 
is introduced to keep the area constant. 



6.4 Parameters 

In experiments with vesicles in shear flow, inertial effects 
are negligible since the Reynolds number Re = jpR* 2 /rj out 
is very small. We express our results using the reduced 
area A* = A /ttR* 2 , defined in Eq. (@|, and the reduced 
shear rate x = T^out^o/K, see Eq. (I3U1) . as relevant di- 
mensionless quantities. 

We set a = 7r/4, n c — 10, and I — 0.008a. This implies 
a viscosity ?7 out = rji n ~ 28.0^/m s kBT/a. We use L x — 
150a, L y = 90a, R* = 15.3a, and v wa u such that Re < 
0.1 for all the cases we considered with 0.5 < \ < 10.0. 
Finally, we set m p = 10m s , N p = 96, At p = At s / '20, r p = 
0.9a, r = a, k = A0k B Ta, k A = Q.5k B T, k h = 4000fc B T. 
The area Ao is chosen in such a way that 0.7 < A* < 0.95. 
With the choices for k A and kh the area and the length 
of the vesicle are kept constant with a deviation of less 
than 1% of the target values for all simulated systems. A 
snapshot of a simulated vesicle and the resulting velocity 
field for the reduced area A* — 0.95 and reduced shear 
rate \ = 5.6 is shown in Fig. [TJ 



7 Results and discussion 



V,; - 2- 



m s + m 



p 



Vj + 2 - ^ (vi-v,-)- 



m s + m 



(103) 



p 



To avoid that a solvent particle moves too far inside a 
disk, we require that / <C r p . The collision step (|98| is 
performed only on those solvent particles which did not 
scatter. If the collision step were executed also on the scat- 
tered solvent particles, they might continue to collide with 
the same disk in the next time step. The fluids in the in- 
terior and exterior of the vesicle are taken to be the same, 
in particular to have the same viscosity. 

A chain of disks of finite radius r p has an inner length 
available to the solvent particles which is smaller than the 
outer length. Since the solvent has the same density inside 
and outside, the outer fluid exerts a compression force on 



In Fig. [2] we show the stationary deformation correlations 
(<5a^), (5b 2 n ) as a function of the mode number m for 
A = 0.163 and \ — 9-3. We also show a fit of these cor- 
relations for m > 3 with the theoretical prediction ([65]) . 
From the fit we can extract the tension S. In this partic- 
ular example we obtain Z^t — 103k/-Rq, whereas theory 
predicts Zthcor - 113/c/i?§ from Eqs. (f6B]> - (f68|) . 

The mean field treatment (|6"5")) predicts (<5a^) = (<5o^ n ) 
for all m. We can see that this holds only for m > 3. As 
explained in Sec. 14. 3[ this is due to the fluctuations in the 



line tension. In the inset of Fig.[5J we compare (Sa 2 ,), (5b 2 



with the low temperature expansion 
agreement. 



2/5 \ uu 2l 

with very good 



7.2 Tension vs shear rate 

Fig.[3]shows the extracted dimcnsionless tensions EmR 2 / '« 
for different dimcnsionless shear rates < x < 10 ancl f° r 
two different excess lengths A = 0.163 and A ~ 0.340. 
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The agreement with the theoretical prediction from a nu- 
merical solution of Eqs. (|66|) -([68 |) is satisfactory. The fact 
that this function is nearly a straight line implies that 
the large-shear-rate approximation (|14p is valid down to 
small shear rates. We find a crossover shear rate Xc, below 
which there are deviations from a linear behavior. Theo- 
retically, Eq. (|7Tj) gives an order-of-magnitude estimate of 
Xc ^ 3.09 for A = 0.163 and X c ^ 0.99 for A = 0.340. 

7.3 Autocorrelation function 

In Fig. HJ the time autocorrelation function (03(^)03(0)) is 
shown as a function of dimensionless time r yt. The data fol- 
lows the expected exponential decay (|64[) very well. From 
the amplitude (5a^(0)5a^(0)} we can extract a tension 
S ~ 65k/ Rq, while from the time constant we deduce 
£ ~ 44K/i?Q. Theory predicts -Stheor ~ 73k/Rq. Given 
the rather noisy data, this agreement seems reasonable. 
For moderate shear rates the autocorrelation function has 
decayed before the oscillations implied by Eq. (|64|) become 
noticeable. Even for the large shear rate x — 5.6 used in 
Fig-Hl the oscillations are barely visible. For the same rea- 
son the build up of cross correlation (a(t)b(0)) is hidden 
in the numerical noise. 



7.4 Inclination angle 

We compare the averaged inclination angle (0) for dif- 
ferent reduced areas A* with Eq. (|78|) in Fig. El valid in 
the quasi-circular limit. The agreement with simulation 
data is satisfactory, given the large error bars. The 2d 
Keller-Skalak theory, which is also shown in the plot, gives 
slightly better agreement. 

In Fig. [6j we show the fluctuations of the inclination 
angle (AO 2 ) = (0 2 ) — (0) 2 as a function of the shear 
rate. The theoretical scaling is given by Eq. (|86| . and is in 
excellent agreement for A* = 0.85 (A ~ 0.53), A* = 0.9 
(A ~ 0.34), and A* = 0.95 (A ~ 0.163). The scaling of the 
fluctuations of for A* differs significantly for A* = 0.7 
(A ~ 1.23). In the deterministic vesicle with such 

a low reduced area would tumble within the quasi-circular 
theory. This implies that the quasi-circular approximation 
works well for A < 0.53 (corresponding to A* > 0.85) in 
two dimensions. 



7.5 Tank-treading frequency 

Finally, we show the rescaled tank-treading frequency uitt/j 
as a function of A* in Fig. [7] Again, agreement with the 2d 
Keller-Skalak theory is quite good, while the quasi-circular 
theory neglects the effect of the vesicle shape on the flow 
and would predict w t ^ C /j = 1/2, see Eq. ([79"1) . 

8 Summary 

We have studied the fluctuations and deformation of a 2d 
vesicle in shear flow at finite temperature. In the limit of 



Fig. 1. Snapshot of the vesicle and the velocity field around 
it, taken from simulation data for reduced area A* = 0.95 and 
reduced shear rate \ — 5-6 (see Eq. (|39[l '). The disks represent 
the beads forming the membrane and are plotted to scale. 

small deformations from a circle, we have derived analyt- 
ical Langevin-type equations of motion, which are nonlin- 
ear due to the length constraint. A mean-field treatment 
allows approximate predictions for the stationary correla- 
tion functions and time autocorrelation functions of the 
deformation amplitudes, which agree quantitatively with 
simulation data. Deviations of the stationary correlations 
from the mean-field predictions in the lowest mode are 
explained quantitatively in a low temperature expansion 
of the original constrained Langevin equations. The mean 
inclination angle and the tank-treading frequency are bet- 
ter described by a deterministic 2d Keller-Skalak theory. 
Fluctuations of the inclination angle are also determined 
quantitatively. Theory and simulations agree well for low 
excess lengths, but differ for larger excess lengths. 

The good quantitative agreement of mesoscale simula- 
tions of vesicles in flow with detailed theoretical calcula- 
tions demonstrates the predictive power of these simula- 
tion methods for more complex flow geometries. 
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